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MULTI-DIMENSIONAL  SIGNAL  PROCESSING  RESEARCH  PROGRAM 


1.  INTRODUCTION  AND  SUMMARY 

The  Lincoln  Laboratory  Multi-Dimensional  Signal  Processing  Research 
Program  was  initiated  in  FY  80  as  a  research  effort  directed  toward  the 
development  and  understanding  of  the  theory  of  digital  processing  of  multi¬ 
dimensional  signals  and  its  applications  to  real-time  image  processing  and 
analysis.  A  specific  long-range  application  is  the  automated  processing  of 
aerial  reconnaissance  imagery.  Current  research  projects  which  support  this 
long-range  goal  are  image  modeling  for  segmentation,  classification,  and 
target  detection,  techniques  for  adaptive  contrast  enhancement,  and  multi¬ 
processor  architectures  for  implementing  image-processing  algorittaa.  In 
this  Semiannual  Technical  Summary  we  shall  describe  in  detail  our  work  to 
date  on  the  detection  of  anomalous  areas  using  an  image-modeling  approach. 

In  the  area  of  adaptive  contrast  enhancement  we  have  recently  completed 
the  initial  development  of  a  software  module  for  implementing  the  algorithms 
previously  developed  under  this  program.  This  software  module  was  written  in 
the  programming  language  "C"  to  facilitate  its  transfer  to  the  RADC 
Automatic  Feature  Extraction  System  (AFES),  which  utilises  a  UNIX-based 
operating  system.  We  are  currently  in  the  process  of  transferring  this 
software  module  and  its  documentation  to  RAUC/XRRE  for  testing. 

We  have  begun  the  paper  development  of  a  etrawaan  mult i-procesaor 
architecture  for  image-processing  applications.  We  have  been  looking  at 
generic  image-processing  operations  such  as  linear  filtering,  median 
filtering,  nonlinearities,  segmentation,  edge  detection,  and  Fourier 
transformation.  Based  on  operations  such  as  these,  it  seems  that 
straightforward  data  partitioning  techniques  can  be  used  in  most  cases  to 
fit  the  operation  into  a  mult i-processor  implementation.  We  plan  to 
continue  refining  our  architectural  approach  during  the  remainder  of  FY  82  to 
incorporate  other  image  operations  and  to  address  the  problems  of 
programmability ,  I/O,  inter-processor  r ^-mirations,  and  required 
performance  levels. 
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The  remainder  of  this  report  is  devoted  to  a  comprehensive  treatment  of 
both  our  theoretical  and  empirical  results  in  the  area  of  target  detection 
within  aerial  reconnaissance  imagery. 

2.  TARGET  DETECTION  BY  PARAMETER  TRACKING 

This  component  of  our  research  relates  to  the  problem  of  detecting 
targets  (i.e.,  anomalous  areas)  in  aerial  photographs  and  is  an  outgrowth  of 
our  work  in  image  modeling,  segmentation,  and  classification.  We  can  loosely 
define  the  target-detection  problem  as  the  detection  of  man-made  objects  in  a 
textured  background  (e.g.,  trees,  grass,  fields,  etc.). 

Usually,  in  detection  theory  the  target  (or  signal)  is  added  to  the 
background  (or  noise),  and  filtering  procedures  are  well  established  for 
increasing  the  aignal-to-noiae  ratio.  In  image  processing,  however,  the 
target  pixels  replace  the  background  pixels.  Motivated  by  this  observation, 
we  model  an  aerial  photograph  as  s  2-D  random  process  which  is  the  output  of 
a  continuously  space-varying  linear  filter.  We  assume  that  a  textured 
background  (i.e.,  the  2-D  filter  representing  the  background)  is  slowly 
varying,  while  an  anomalous  area  of  the  image  reflects  a  quick  change  in  the 
filter  characteristic*. 

One  possible  interpretation  of  the  target-detection  problem  then 
consists  of  detecting  large  changes  in  2-D  filter  parameters.  That  is,  we 
"track14  the  parameters  of  the  model;  when  a  target  area  is  reached,  the 
parameters  tend  to  change  quickly,  indicat ing  the  presence  of  a  target. 
Furthermore,  we  shall  assume  no  knowledge  of  the  target  or  background  except 
that  the  target's  spectral  characteristics  differ  from  that  of  its 
background.  In  particular,  the  target  say  have  high  or  low  frequencies 
compared  with  the  background  and/or  a  different  mean  or  variance  than  the 
background.  This  assumption  differs  from  that  of  many  available  target- 
detection  algorithms  which  rely  on  only  a  change  of  intensity  and/or 
variance  of  the  target  from  the  background. 
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Our  preliminary  results  in  this  approach  to  target  detection  are 
discussed  in  the  previous  Semiannual  Technical  Summary,1  where  we  addressed 
this  problem  in  one  dimension  to  develop  the  theory  and  our  intuition  before 
extending  the  results  to  the  two-dimensional  problem. 

In  Sec.  2.1,  we  shall  review  and  elaborate  on  this  one-dimensional 
theory.  In  particular,  we  derive  nonrecursive  and  recursive  least-squares 
procedures  for  parameter  tracking  and  detection  when  the  data  are  seen 
through  a  finite-extent  sliding  window.  From  these  derivations,  we  see  that 
the  recursive  approach  is  preferable  (computationally)  and  also  gain  further 
insight  into  the  results  of  Ref.  1.  Furthermore,  these  results  lead  us  to 
some  extensions  of  Ref.  1  for  the  one-dimensional  problem  which  appear  to 
hold  promise  for  the  two-dimensional  case.  Specifically,  we  introduce  the 
use  of  prediction  error  in  detecting  anomalous  areas  of  a  signal.  A  target 
is  detected  when  a  large  error  occurs  in  predicting  a  sample  of  the  signal 
from  its  past  (i.e.,  causal  prediction).  Almost  always,  the  prediction  error 
compares  favorably  with  the  change  in  filter  parameters  when  used  in 
detecting  small  anomalous  areas  of  a  signal.  We  also  investigate  the 
prediction  of  a  signal  sample  from  its  surrounding  values  (i.e.,  noncauaal 
prediction).  Such  an  approach  (noncauaal)  to  detection  may  be  particularly 
useful  in  a  2-D  context  where  no  directionality  is  assumed. 

In  Sec.  2.2,  we  extend  a  subset  of  our  l-D  results  to  develop  a  2-0 
target-detection  algorithm  which  relies  on  tracking  Che  parameters  of  a 
continuously  apace-varying  2-D  filter.  As  in  the  l-D  case,  we  explore  both 
recursive  and  nonrecursive  structures  for  least-squares  parameter  estimation. 
Unlike  the  l-D  case,  however,  we  shall  see  that  it  may  be  desirable 
(computationally)  to  perform  2-D  least-squares  directly  (i.e., 
nonrecuraively)  at  each  pixel  rather  than  recursively,  nevertheless,  our 
derivation  of  a  2-0  recursion  lends  insight  into  the  nature  of  filter 
parameter  changes  across  an  image. 

As  before,  detection  of  an  anomalous  area  of  an  image  is  based  on 
changes  in  the  ectimated  filter  parameters.  We  shall  show  that,  because  of  a 
lack  of  image  directionality,  the  choice  of  the  mask  of  the  2-D  filter  model 
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can  play  an  important  role  in  accurate  estimation  of  the  parameters  (and  thus 
their  change  used  in  detection).  We  shall  illustrate  such  variations  on  a 
baseline  algorithm  by  detecting  small  targets  (e.g.,  4x4,  3x3,  and  2x2 
pixels)  within  RADC  image  data.  Finally,  we  comment  on  the  extension  to  2-D 
of  the  use  of  causal  and  noncausal  prediction  error  in  detection.  In 
particular,  we  conjecture  that  the  2-D  noncausal  prediction  error  may 
provide  a  significant  improvement  over  the  present  algorithm  (based  on 
parameter  changes)  in  detecting  small  objects  in  a  cluttered  background  such 
as  a  forested  area.  Consequently,  investigation  of  various  2-D  prediction- 
error  functions  will  be  a  large  component  of  our  future  research  in  target 
detection, 

2.1  Target  Detection  in  One  Dimension 

Here,  we  present  an  algorithm  for  target  detection  in  1-D,  Our 
analysis  is  used  as  a  stepping  stone  to  the  development  of  an  analogous  2-D 
algorithm  in  the  following  section. 

2.1.1  Nonrecursive  Parameter  Estimation 

Suppose  that  a  1-D  sequence  s(n)  follows  an  autoregressive  model,  and 
thus  is  generated  by  a  difference  equation  of  the  form 


s(o)  “  l  a(k)  s(n  -  k)  ♦  w(n)  (1) 

k-1 

where  w(n)  is  xero-taean  white  noise.  Our  objective  is  to  estimate  from  a(n) 
the  model  parameters  a(k)  for  k  e  Jl,P).  Further,  suppose  that  we  have 
available  a  segment  of  s(n)  for  n  in  the  interval  l-P  ♦  ng.njl,  that  ia, 

-  no  *  P  ♦  1  data  points.  We  then  define  the  error  e(n)  over  the 
interval  (ng.nj)  (which  we  assume  is  N  points  in  duration,  i.e.,  N  » 
nj  -  oq  ♦  1)  by 


e(a)  ■  s(n) 


a(k)  s(n  -  k) 


n  e 


k-1 


(Vni! 


(2) 
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e(n)  can  be  interpreted  as  the  error  in  predicting  s(n)  from  its  past 
P  samples.  Our  goal  then  becomes  to  minimize  the  sum  of  the  squared 
prediction  errors  given  by 


To  accomplish  this  minimization,  we  first  place  the  problem  in  a  matrix 
algebra  framework.  In  particular,  we  define  the  vectors  0  and  a  by 


and  the  matrix  S  by 


Thus,  the  error  Etnj]  in  Eq.  (3)  becomes 

Efftjl  ■  [Sa  -  o)T  tSa  -  a}  (4) 
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where  T  denotes  matrix  transpose.  Differentiating  Eq.  (A)  with  respect  to 
a  to  find  the  least-squares  error,  we  obtain2 

a  -  (STS]  1  STo  (5) 

T 

where  S  S  is  a  P  x  P  matrix  (P  being  the  order  of  our  model) .  The  solution 

Eq.  (5)  is  often  referred  to  as  the  covariance  method  of  linear  prediction. 

Thus,  the  inverse  of  a  P  x  P  matrix  is  required  to  compute  a.  This  will 

require  on  the  order  of  P3  operations.2 

Suppose  now  that  we  wish  to  obtain  an  estimate  of  a  for  many  N-point 

segments  and,  in  particular,  over  the  intervals  [n0  ♦  m,  n1  +  m)  where  m  ■  0, 

1,2,...  .  In  the  next  section,  we  derive  a  method  of  updating  the 

T  -1 

coefficients  a  without  recomputing  the  matrix  (S  S) 

2.1.2  Recursive  Parameter  Estimation* 

Let  us  suppose  that  a  new  data  point  is  given  (i.e.,  s(nj  +1)]  and  that 
one  is  eliminated  (i.e.,  s(no)1  an  that  our  new  data  set  runs  from  ng  -  P  ♦  1 
to  ni  ♦  1.  The  sum  of  the  squared  errors  over  the  interval  (n0  *  1,  n^  ♦  1] 
then  becomes 

Elnj  ♦  1]  «  Elt^)  ♦  e2(nj  ♦  1)  -  e2(nQ)  (6) 

where  ♦  1)  is  the  error  in  predicting  the  new  sample  sin^  ♦  i): 

P 

e(«l  ♦  1)  »  s(n,  >  1)  -  I  a(k)  s(n,  ♦  1  -  k)  (7a) 

1  k-1  1 


*The  derivation  in  this  section  parallels  that  of  Ref.  3,  where  the  data  act 
increases  with  time. 
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and  e(ng)  is  Che  error  in  predicting  the  old  sample  s(ac): 


e(n  )  -  8(nn)  -  £  a(k)  s(nn  -  k) 

U  k-1 


In  matrix  form,  we  can  write  Eq.  (6)  as 


where 


E[n1  ♦  1]  -  E[nj]  ♦  [da  -  s(nj  ♦  l)}2  -  [ca  -  s(n0)l2 


s(nt  ) 
s(ni  -  1) 


a(ni  ♦  1  -  P). 


sCuq  -  1) 
s (no  -  2) 


a(no  *  P) 


Differentiating  Eq.  (8)  with  respect  to  a,  we  obtain2 


a  *  [S^S  -  cTc  ♦  d^dj  ^  (S*^0  *  d^sCn^  +  1)  -  c^»(o^)J 


How,  we  can  simplify  Eq.  (9)  as 


where 


a  •  tt  ♦  cVl  (STo  ♦  GTQ] 


sCn^  +  1) 


-s(n0) 


and  where  R  *  S  S.  We  wish  to  avoid  direct  computation  of  the  matrix  inverse 
in  Eq.  (10).  To  accomplish  this,  we  rely  on  the  matrix  identity  given  by 

[A  +  BC]-1  »  A-1  -  A~1B  [CA-1B  +  I]-1  CA_1  (11) 

T 

where  I  is  the  identity  matrix.  Letting  A  =  R,  3  *  G  ,  and  C  =  F  in 
Eq.  (10)  ,  we  obtain 

a  =  [R-1  -  R~1GT(FR"1GT  +  I)"1  Fr"1]  [STo  +  GTQ]  .  (12) 

Note  that  R“*  is  assumed  known  from  the  previous  computation  of  a*  At  this 
point,  for  clarity,  let  us  change  notation  so  that  a  becomes  a[ni  ♦  1], 
i.e.,  the  filter  parameters  estimated  within  the  interval  [n0  ♦  1,  nj  +  1], 
Then,  Eq.  (12)  can  be  written  as 

-IT  -IT  -IT  -1  ~1  T 

alnj  +  1]  -  R  S  o  -  R  G  (FR  G  +1)  FR  So 

-IT  -IT  -IT  -1  -1  T 

+  R  GQ-R  G(FR  G  +  I)  FR  G  Q  .  (13) 

-IT 

Noting  that  a[n^]  ■  R  So  and  defining  a  "gain  matrix"  K  as 
-IT  -IT  -1 

K  “  R  G  (FR  G  +1)  (14) 


we  have 
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a[n1  +  1]  -  afnj]  -  XFafn^  +  R_1GTQ  -  KFlfVq 

-  aCn^  -  KFa[nj3 

+  R"1GT(FS”1G1  +  I)-1  (FR_iGr  +  I)Q  -  KFr"1GTQ 
-IT  -IT 

-  alnj]  -  KFa[n1)  +  K(FR  G  +  l)Q  -  KFR  G  Q 

-  alnj]  -  RFafn^  +  KQ  .  (15) 


Finally,  we  have 

^ [n^  +  1]  »  a[n^]  +  K(Q  -  Fa[n^])  (16) 

which  provides  a  recursive  means  of  updating  a[n]. 

Now  let  us  consider  the  gain  matrix  K.  The  matrix  computation  involved 
in  finding  K  can  be  represented  as 


(17) 


Thus,  we  need  to  compute  the  inverse  of  a  matrix  (i.e,,  FR  G  +1)  which  i« 

of  else  2x2,  in  contrast  to  the  direct  approach  which  requires  the  inverse 

,  t  -1 
of  a  P  x  P  matrix  [entailing  0(PJ)  operations].  Note  from  Eq.  (12)  that  R 

can  also  be  recursively  updated  as 


-  KFR 


(18) 


which  requires  some  additional  matrix  multiplications  of  order  P2 . 
Consequently,  the  bulk  of  the  computation  Ilea  not  in  the  matrix  inversion, 
but  in  matrix  multiplications  of  order  P2. 
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2.1.3  Detection  Functionals 

It  was  pointed  out  in  the  introduction  to  Sec.  2  that  one  procedure  for 
detecting  anomalous  areas  of  a  signal  is  to  detect  abrupt  changes  in  the 
parameter  estimates  of  the  time-varying  filter.  If  the  model  parameters  of  a 
stationary  background  are  known,  then  one  obvious  choice  of  a  detection 
functional  is  the  mean  squared  difference  between  the  known  background 
parameters  and  the  parameter  estimates.  If  the  data  deviate  from  the 
background  (i.e.  ,  a  target  is  present),  this  mean  squared  error  should 
suddenly  change  to  account  for  the  large  error  that  would  be  encountered  in 
predicting  the  first  few  points  of  the  new  data  set  if  a[ni  +  1]  didn't 
change  much  from  a[nj].  Examples  of  such  abrupt  changes  in  this  functional 
were  illustrated  in  kef.  1. 

This  section  describes  and  demonstrates  with  examples  alternative 
detection  functionals  which  do  not  require  background  information  or 
assumptions  of  stationarity. 


Thr  Coefficient  Change  Functional:-  To  free  ourselves  from  a  priori 
information  about  the  background  and  the  atationarity  assumption,  we  first 
considered  a  detection  functional  of  the  form: 


C[nL] 


lUj)) 


♦  1]  -  aUj]) 


which  is  the  mean  squared  difference  between  two  consecutive  parameter 

estimates,  fe’e  might  expect  C[nil  to  jump  only  when  our  sliding  window  first 

hits  the  target  and  when  the  window  just  leaves  the  target.  Something 

similar  to  CfniJ  was  considered  by  Dove  and  Qppenhcim4  who  referred  to 

Eq.  (19)  as  the  coefficient  change  functional. 

Although  this  functional  has  been  used  with  good  success  in  both  1-D  and 

in  extensions  to  2-D,  a  potential  problem  exists.  This  change  functional  is 

sometimes  not  responsive  to  small  extent  targets  in  a  highly  cluttered 

background  (e.g. ,  trees).  Note  that  the  parameter  estimates  rely  on  the 

T 

correlation  function  R  ‘  S  S.  A  small  extent  target  may  not  change  R 


sufficiently  to  generate  a  corresponding  change  in  filter  parameters. 
However,  such  absence  of  change  in  the  parameters  may  actually  aid  a 
detection  functional  based  on  prediction  error  -  as  we  shew  in  the  following 
section. 

The  Prediction  Error  Functional:-  Let  us  return  to  Eq.  (16)  which  we 
new  write  as 

a(n^  +  1]  *  afnj]  +  ketn^]  (20) 


where 


etn^  *  Q  -  Fa[n1  ] 


m 

) 

a> 

S' 

►— 

+ 

_ 1 

d  ‘ 

-s(nQ) 

-c 

s(nL  ♦  1)  -  da(nl 1 
s(n0)  ♦  caln^ ] 

eiKi] 

i 


(21) 


Frcm  Eqa.  (7)  and  (6),  we  see  that  ej(n)  in  Eq.  (21)  ia  the  error  in 
predicting  the  new  signal  ssrple  with  the  old  parameters.  When  this  new 
signal  sample  s(nj  ♦  1)  is  a  member  of  an  anomalous  area,  we  expect  Che  error 
in  predicting  s(’*i  ♦  1)  to  change  abruptly  when  the  old  set  of  coetiicienta 
atn^J  is  baaed  on  data  not  encowpacaing  the  target.  (tote  that,  although 

Che  coefficient  change  a[ui  ♦  l]  -  a(ni)  is  a  function  of  ei(n)  (i.e., 
a(nj  ♦  1)  -  ajnj]  ■  Kelnjj),  the  time-varying  gain  term  K  and  e2(o)  may  drown 

out  the  change  in  ei(n). 
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It  is  possible  to  generalize  this  notion  of  prediction  error  for  target 
detection.  In  particular,  e^Cn)  in  Eq.  (21)  represents  a  causal  prediction 

error,  i.e. ,  s(ni  +  1)  is  predicted  from  the  past  P  values  of  s(n).  Clearly, 
it  is  possible  to  formulate  the  least-squares  estimation  problems  of 
Secs.  2.1.1  and  2.1.2  in  terms  of  noncausal  prediction  of  a  sample  of  s(n) 
from  its  surrounding  neighbors,  i.e.,  a  smoothing  of  s(n)  from  points  to  its 
left  and  right: 


8(n)  *  £  a(k)  s(n  -  k)  +  w(n) 

&oQ 


(22) 


We  have  applied  such  noncausal  detection  to  1-D  with  success.  Noncausal 
prediction  may  be  particularly  useful  in  2-D  detection  where  there  exists  no 
apparent  directionality.  In  the  next  section,  we  compare  by  example  the 
various  detection  functionals  discussed  in  this  section. 

2.1.4  Examples 

Consider  a  sequence  s(n)  of  the  form: 

s(n)  -  0.95s(n  -  1)  ♦  w(n)  (23) 

where  w(n)  is  zero-mean  white  noise.  A  sample  function  of  s(n)  is  shown  in 

Fig.  1.  A  4-point  high-frequency  target  at  n  *  90,  91,  92,  and  93  is  also 

shown  there.  Figure  2(a)  depicts  the  coefficient  estimate  based  on  a 

16-point  vindw,  and  Fig.  2(b)  shows  the  corresponding  coefficient  change 

C[n).  Note  that  a  global  change  in  mean  and  other  local  phenomena  in  the 

background  has  been  emphasized  over  the  target.  Figure  2(c)  depicts  the 

2 

squared  prediction  error,  [n ] ;  the  target  is  clearly  detected. 

In  a  second  example,  we  analyzed  a  cross  section  of  pure  tree  data  from 

the  field-tree  RADC  image.  Figure  3  depicts  this  cross  Section  with  its  mean 

removed  (see  Sec. 2. 1.5),  with  a  3-point  constant-level  target  imposed  at  n  * 

64,  65,  and  66.  A  5-pole  model  was  assumed  and  a  16-point  sliding  window 

2 

was  used.  C(n)  and  e. (n)  are  shown  in  Figs.  4(a)  and  (b),  respectively. 

2  * 

Note  that  e^tn)  is  without  the  false  alarm  of  C(n). 
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Consider  again  the  sequence  generated  by  «5q.  (23)  in  Fig.  it  but  now 
where  we  assume  a  model  of  the  fora 


s(n)  •  s(l)  s(n  -  1)  ♦  s(2)  s(n  ♦  1)  ♦  w(n) 


(24) 


i.e,,  s(n)  is  s  function  of  its  two  nearest  adjacent  neighbors.  Figure  S 
depicts  the  noncausal  prediction  error  corresponding  to  the  signal  in  Fig.  1 
and  which  uses  a  16-point  sliding  window.  Finally,  the  noncausal  model 
Eq.  (24)  was  used  to  generate  a  prediction  error  from  the  image  cross  section 
given  in  Fig,  3.  The  prediction  error  is  depicted  in  Fig.  6.  The  target  is 
clearly  detected  in  both  cases. 

2.1.5  Removal  of  the  Mean  Level 

It  should  be  noted  that  the  image  cross  section  of  Fig.  3  contains  s 
noncero  mean  which  violates  our  model  Eqs.  (1)  and  (22)  where  seronseaa  white 
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e.  In]  (noncaua*!)  for  Fig.  3. 


noise  is  assuned  to  drive  a  linear  system.  One  approach  to  addressing  a 
time-varying  nonzero-mean  level  is  to  incorporate  it  within  our  model  so 
that  we  estimate  the  mean  as  well  as  the  parameters  a(k): 


P 

sCn)  “  £  a(k)  8(n  -  k)  +  w(n)  +  C  (25) 

k-1 


where  C  is  a  mean  level  added  to  the  input  w(n).  The  output  will  then 
contain  a  nonzero  mean. 

Now,  if  we  perform  least-squares  estimation,  we  find  that  the  best 
estimate  of  the  a(k)'s  is  given  approximately  by  (this  derivation  parallels 
that  in  Sec.  2.1.1) : 


a 


/  T  2.-1  .  T  2. 
(S  So  -  i  )  (S  o  -  m  ) 

S  8 


where  N  is  the  data  length  and 


(26) 


m 

s 


"1 

I 

n-n 


0 


a(n)/N 


(27) 


which  is  an  estimate  of  the  mean  of  the  output.  The  best  estimate  of  C  is 
given  approximately  by 


C  *  m  1 

a 


P 

l  a(k) 


k-1 


(28) 


Equation  (26)  simply  implies  that  we  first  estimate  the  mean  of  a(n)  through 
Eq.  (27),  subtract  it  from  s(n),  and  then  solve  Eq.  (5)  of  Sec.  2.1.1.  Tb 
find  C,  we  "inverse  filter"  the  mean  at  the  output  to  obtain  the  mean  at  the 
input. 

Removal  of  the  mean  may  not  only  lead  to  more  accurate  determination  of 
the  a(k)*s,  but  also  provide  an  Additional  parameter  in  our  functional  for 
target  detection.  We  shall  discuss  this  possibility  in  the  next  section. 
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2.2  Target  Detection  in  Two  Dimensions 

One  obvious  approach  to  extending  our  1-D  results  to  2-D  target 
detection  i9  to  apply  1-D  recursive  least-squares  estimation  along  rows  or 
columns  of  an  image.  Of  course,  such  a  procedure  is  not  "optiIQal,,  since  the 
correlation  of  a  sequence  s(n,m)  is  considered  in  only  one  direction.  We 
have  found  empirically  two  problems  with  this  approach:  (1)  various 
detection  functionals  used  have  a  large  dynamic  range  (i.e.,  the  functional 
exhibits  little  consistency  from  line  to  line),  and  (2)  an  undesirably  long 
window  length  (along  each  row  or  column)  may  be  required,  resulting  in 
excessive  Hraamory"  of  the  functional.5'5 

In  this  section,  we  describe  a  least-squares  procedure  which  invokes 
rows  and  columns  simultaneously,  i.e.,  we  address  the  2-D  problem  directly. 
After  considering  both  nonrecursive  and  recursive  estimation,  we  devise  a  2-D 
target  detection  algorithm  which  is  a  significant  improvement  over  its  line- 
by-line  counterpart.  We  demonstrate  the  algorithm  by  applying  it  to  the 
detection  of  anomalous  areas  within  aerial  photographs  obtained  from  RADC. 

2.2.1  Nonrecursive  Parameter  Estimation 

Conceptually,  2-D  least-squares  estimation  is  similar  to  the  1-D  case. 

We  begin  with  a  2-D  sequence  which  follows  an  autoregressive  model,  so  it  is 
given  by  a  difference  equation  of  the  form 

s(n,m)  •  l  £  *(j,k)  s(n  -  j,  a  »  k)  ♦  v(n,o)  (29) 

i  R 

(j,S;)*(0,0) 

where  v(n,o)  is  eero-mean  white  noise.  For  the  moment,  we  shall  assume  that 
the  prediction  coefficients  a(j,k)  fall  within  s  (P  x  Q)  first-quadrant 
mask,  as  illustrated  in  Fig.  7.  For  simplicity,  we  limit  our  derivations  to 
this  class  of  prediction  masks,  although  it  is  clearly  applicable  to  more 
general  mask  shapes.  In  particular,  we  shall  investigate  second-quadrant  and 
nonsymmetric  half-plane  masks  in  Sec.  2.2.5. 

Our  objective  is  to  estimate  from  s(n,m)  the  model  parameters  a(j,k)  for 
j  *0,1,...,?  and  k  ■  0,1,..., Q,  with  j  ■  k  *  0.  Further,  let  ua  suppose 


19 


that  we  have  available  a  square  block  of  data,  i.e.  ,  s(n,m)  for  (n,ra)  e  [-P  ♦ 
ng.n^]  x  { — Q  ♦  mg  .idj  ]  (see  Fig.  8).  We  then  define  the  error  e(n,m)  over  the 
region  I,  given  by  I  *  [no.nj.]  x  (mo, mi],  uQ 

p  Q 

e(n,m)  ■  s(n,m)  -  £  £  a(j,k)  s(n  -  j,  m  -  k)  (n,m)  e  I  .  (30) 

j«0  k»0 
j-k*0 

Our  goal  becomes  to  minimize  the  sum  of  the  squared  errors  given  by 


2 

E(n,,m,]  «  [  Is  (n,m)  .  (31) 

n-n0 

The  approach  we  take  ia  to  transform  the  2-D  problem  to  a  1-0  problem  so 
that  a  l-D  least-aquarea  solution  ia  applicable.  Note,  however,  that  we  will 
still  have  solved  the  2-0  least-squares  problem.  In  particular,  we  wish  ec 
transform  Eq.  (31)  into  a  form  similar  to  the  1-0  error  expression  in 
Eq.  (/>).  To  accomplish  this  transformation,  we  define  the  following  vectors 
a(nj  ,ra ^ ]  and  o  by 
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a( u^  ,m^ 


'  a(0,l)" 

4 

8(n0,m0) 

it 

a(0,2) 

8(n0,  mQ  +  1) 

* 

• 

• 

• 

• 

• 

a(0,Q> 

sCnQ.n^) 

a(l ,0) 

s(n0  +  1,  m0) 

a(l ,1) 

8(pq  +  1,  mQ  +  i) 

• 

• 

• 

c 

(i 

3Q  -  1)  0  “ 

• 

(P 

a(l,Q) 

8(nQ  +  1,  n^) 

I 

. 

• 

• 

♦ 

• 

• 

a(P>0) 

a(npfflg) 

a(P,l) 

s< n j, ,  mQ  <■  1) 

• 

• 

9 

• 

• 

• 

1 - 

to 

O 

( 

a(a1,ia1) 

and  the  matrix  S  by 


where 


—  FQ  '  1 

(•(n,;  » 

j  -  n,  «0  -  r-...«(n0  *  j  -  n,  «„  ♦  o)l 

...  lila. 

(*<"o  ♦ 

j  -  0,  «0  ♦  •  -  1)...»(q0  ♦  i  -  0,  «0  ♦  i  ♦  0)1 

...  (•<«„ 

i'<n0  ♦ 

j  -  O,  ■l  -  0)1 

...  I«(n0 

(33b) 


and  where  we  have  assumed  the  known  data  segment  to  be  of  extent  N  x  N.  Note 
that  o  is  a  vector  consisting  of  the  concatenation  of  the  rows  of  s(u,m)  over 
I,  afn^.m^]  ia  a  vector  consisting  of  the  concatenation  of  the  rows  of 
a(j,k)  for  (j,k)  e  [0,P]  x  [0,Q]  with  j  ■  k  *  0,  and  £  is  a  matrix  which 
consists  of  the  concatenation  of  rows  of  varioua  sub-sequencea  of  the  known 
s(n,uk)  required  in  predicting  each  value  of  s(n,m)  over  I. 

Therefore,  we  can  write  Eq.  (31)  as 


n  1  "I 


Eln^mJ  -  l  l  e^n.o) 


n“Q0  °“a0 


T 

(Safttj  ,n^)  -  o)  (Sainj.o^]  -  cr) 


From  the  similarity  of  Eq.  (34)  with  Eq,  (4),  we 
minimising  Eq.  (34)  with  respect  to  s(nj,mi]  as 


-I  T 

ajn^.o^)  ■  R  So 


(34) 

can  write  the  solution  to 


(35a) 


where 


T 

a  -  s  s 


(35b) 


Note  that  the  matrix  R  is  of  extent  (PQ  -  1)  x  (PQ  -  1)  and  difficulty  in  its 
inversion  is  dependent  on  the  model  order,  not  on  the  size  of  the  known  block 
of  data. 

In  particular,  since  R  is  generally  not  Toeplts,  its  inversion  will 

require  on  the  order  of  (PQ)3  operations  (i.e.,  multiplications  within,  for 

example,  Cholesky's  decomposition  algorithm2).  It  should  be  pointed  out  that 

the  computation  of  R"1  can  probably  be  reduced  by  considering  its  block-like 

structure  resulting  from  the  conversion  of  a  2-D  problem  to  a  1-D  problem. 

Thus,  assuming  P,Q  <<  N,  the  bulk  of  the  computation  is  embedded  within 
T 

forming  R  *  S  S  which  requires  on  the  order  of  N1*  operations. 

2.2.2  Recursive  Parameter  Estimation 

Let  us  suppose  new  that  we  wish  to  estimate  the  parameters  a[ni,rai]  over 
a  large  data  set.  Our  goal  in  this  section  is  to  devise  a  recursive 
procedure  for  updating  the  parameter  estimates  over  the  2-D  plane  without 
recomputing  the  inverse  of  R.  Unlike  the  1-D  problem,  many  ways  exist  to 
recurse  in  the  2-D  plane,  i.e.,  vertically,  horizontally,  diagonally,  etc. 

In  all  these  procedures,  entire  rows  and/or  columns  of  data  enter  and  exit 
the  view  of  a  2-D  sliding  window.  An  alternative  to  this  class  of  recursions 
is  the  class  of  "snake  recursions"  which  allow  only  one  or  a  few  signal 
sonples  to  enter  and  exit  the  window's  field  of  view.  Such  recursions  may 
not  only  reduce  computation,  but  provide  greater  resolution  in  detecting 
small-extent  2-D  targets.  For  the  present,  however,  we  have  investigated  the 
row  (column)  recursion,  and  its  derivation  is  the  primary  emphasis  of  this 
section.  In  the  conclusion  of  this  section,  however,  we  briefly  describe  one 
point-by-point  recursion,  its  computational  advantages,  and  its  potential  for 
high-resolution  detection. 

Now  let  us  suppose  that  we  add  a  new  row  of  data,  i.e.,  a(ni  ♦  1,  m)  for 
m  e  (~Q  ♦  ®q ,  J  and  eliminate  the  known  row  at  n  ■  -P  ♦  n0  (see  Fig.  9). 

Our  goal  is  to  generate  the  coefficient  estimates  a[ni  ♦  1,  mi)  based  on  the 
new  row  of  data  and  the  old  coefficient  estimates  The  now  error 

of  interest  is  given  by: 
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Fig.  9.  Data  update  used  in  2-D  recursive  least  squares. 


E[n1  +  1,  «  ElnpO^j  -  £  e2(nQ,o)  ♦  £  e2(u.  ♦  1,  n)  (36) 

b-Oq  «b 


where  eCnj  ♦  l,  a)  is  the  error  in  predicting  the  new  row  s^  ♦  1,  m) : 


P  Q 

e(n.  ♦  1,  m)  •  s(n.  +!,»)-  £  l  a(j#k)  s(n  ♦  1  -  j,  a  -  k)  (37) 

j*0  k"0  1 

(j,k)*(0,0) 


and  where  e(ng>m)  is  the  error  in  predicting  the  old  row: 


P  Q 

a(n  ,o)  -  s(n  ,a)  -  l  l  a(j,k)  s(n  -  j,  a  -  k)  .  (38) 

j-0  k-0  u 

( j»k)^(0,0) 
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In  matrix  form,  we  can  write  Eq.  (36)  as 


E[n^  +  1;.  mj  -  Etn^mJ  +  (Dat^  +  1,  mj  -  g)  (Da[n1  +  1,  m^]  -  g) 


-  (Ca[n1  +  1,  n^]  -  f)  (Ca[nj  +  1,  m^]  -  f) 


where  a[n^  +  1,  mj]  is  given  in  Eq.  (32),  and  where  g  and  f  are  the  vector 
representations  of  the  old  and  new  rows,  respectively: 


s(uq ,mQ) 

s(n0,  m0  +  l) 


s(n1  +  1,  m q ) 

8(n1  +  1,  mQ  +  1) 


s(nQ  .n^} 


s(n1  +  1,  m  J 


and  where  C  ■  and  D  ■  with  the  N  x  (PQ  -  1)  matrix  given  in 

Eq.  (33b).  Noting  the  similarity  of  Eq.  (40)  with  Eq.  (8),  we  can  write  the 

solution  +  1,  raj  which  minimises  E[n1  +  1,  mj  as 

a^  +  1,  raj  -  (R  -  CTC  +  D^P1  [STo  +  DTg  -  CTf]  .  (41) 


Now,  we  can  simplify  Eq.  (41)  as 

a(nl  +  1,  raj  -  [R  -  GTFp^  (STo  ♦  GTQ] 


where 


-*PQ  -  1*-  •**PQ  -  1*> 
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and  where 


t 

2N 

I 

Following  our  derivation  in  Sec.  2.1.2,  we  then  have  the  recuraive 
computation  of  a[n^  +  1,  m^]: 

a[n1  +  1,  m^]  ■  atn^ra^]  +  K(Q  -  Faln^m  ])  (43) 

where  K  is  given  by  Eq.  (14)  and  R"*1  can  be  recursively  updated  as  in 
Eq.  (18). 

The  matrix  computation  involved  in  finding  K  can  be  represented  as  (with 
L  -  PQ  -  1): 


-*~L*--*2»**  2N - *• 

K  •  t  [r"’1][gT]|  [pR”1CT  +  ij  2N  (44) 

so  that  K  is  an  L  x  2N  matrix.  The  important  point  here  is  that  the  inverse 

of  a  2N  x  2N  matrix  is  required  to  compute  K,  implying  a  number  of  operations 

on  the  order  of  (N)3  operations  per  imsge  pixel  -  clearly  intolerable  for 

windows  of  large  extent.  Thus,  unlike  in  the  1-D  case,  a  recursion  based  on 

Eq.  (43)  may  entail  more  computation  than  performing  2-D  least-squared 

directly  through  Eq.  (35).  Of  course,  we  hsve  not  considered  sny  special 
-l  T 

structure  in  FR  G  ♦  1  which  asy  help  to  reduce  computstion. 

An  slternstive  to  this  rov-by-row  recursion  (which  we  promised  to 
briefly  describe  before  lesving  this  section)  is  s  "snake-like"  point-by¬ 
point  recursion  illustrated  in  Fig.  10.  To  compute  the  parameter  estimates 
under  the  window  in  Fig,  9,  we  must  perform  N  "little"  recursions  each  of 
which  is  1-D  in  nature,  thus  requiring  the  inverse  of  N  2  x  2  matrices  (see 
Eq.  (14)].  Therefore,  we  are  confronted  with  the  trsde-off  between  one 
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inverse  of  a  large  2N  x  2N  matrix  and  N  inverses  of  little  2x2  matrices  - 
clearly  a  significant  reduction  in  computation.  One  other  potential 
advantage  of  this  point-by-point  recursion  is  that  it  does  not  attempt  to 
predict  an  entire  rear  on  each  step.  An  implication  of  this  property  is  that 
if  target  detection  relies  on  a  prediction  error,  a  point-by-point  recursion 
may  be  more  likely  to  detect  very  small-extent  and/or  closely  spaced  targets 
However,  since  N  recursions  are  required  for  each  pixel,  our  detection 
function  becomes  three-dimensional.  Such  speculation  serves  as  a  starting 
point  for  future  research. 

For  the  present,  however,  we  have  put  aside  all  recursive  solutions  and 
devised  a  detection  algorithm  based  on  direct  computation  of  Eq.  (35). 
Nevertheless,  all  this  discussion  about  recursions  is  useful  since  (as  we 
shall  see  in  the  following  section)  it  has  provided  insight  into  the  nature 
of  changes  in  the  parameter  estimates  across  an  image. 


-\r:* 
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2.2.3  Detection  Functionals 

Clearly,  we  can  extend  our  1-D  detection  fuuctionals  to  2-D.  The 
functional  we  have  chosen  to  use  (for  the  present)  in  the  implementation  of 
2-D  detection  algorithms  is  the  coefficient  change  functional  given  by 

Ctn^]  -  (a(n1  +  1,  -  afn^m^)1  (a[n1  +  1,  -  atn^mj)  .  (45) 

From  Eq.  (43),  CtiipHij]  can  be  written  as 

-  {K(Q  -  Fatn^])}1  (K(Q  -  FaU^])} 


eJnpmJlT  fe^n^rn^ 


e^n^rn^ 


e2  *®i 1 


where  e^n^mj  is  the  error  in  predicting  the  new  row  with  the  old 
coefficient  estimates: 


eitnr“i1  "  -  8 


(47a) 


and  ejlnj.m,]  is  the  cjror  in  predicting  the  old  row  with  the  old 
coefficient  estimates: 


e2lV“iJ  “  f  •  (47b) 

Thus,  C(nj,m^)  resembles  a  weighted  sum  of  these  prediction  errors.  It  is 
important  to  note  that  since  Eqs.  (35)  and  (43)  must  yield  the  same  result 
(based,  of  course,  on  the  same  particular  known  segment  of  s(n,m)],  in  either 
case  (i.e.,  for  either  the  recursive  or  nonrecursive  solution),  CluptUj]  can 
be  written  as  Eq.  (46). 
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Fig.  11.  2-D  detection  algorithm. 
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An  alternative  detection  functional  [analogous  to  the  1-D  function 
ejfnj]  in  Eq.  (21)]  is  the  prediction  error  e^fn^a^].  At  each  pixel  of  an 
image,  we  might  consider  abrupt  changes  in  ei[ni,mi].  If  a  target  happens  to 
fall  along  the  new  row  which  just  enters  the  2-D  (sliding)  window’s  view,  we 
expect  a  large  error  in  prediction  at  the  target  location  along  that  row. 

For  similar  arguments  given  in  Sec.  2.1.3  for  the  1-D  case,  this  prediction 
error  may  be  more  sensitive  to  targets  of  small  extent  than  the  coefficient 
change  function.  Furthermore,  in  the  case  of  images,  noncauaal  prediction 
may  be  more  reasonable  to  attempt  than  the  causal  prediction  given  above. 
Judging  from  our  1-D  experiments  of  Sec.  2.1.4,  detection  based  on  e^[n^,m^] 
(or  its  noncauaal  counterpart)  holds  promise,  and  thus  serves  as  a  possible 
avenue  of  future  research.  However,  in  the  remainder  of  this  report,  we 
shall  investigate  the  use  of  the  coefficient  change  function  for  detection. 

2.2.4  The  2-D  Target-Detection  Algorithm 
Our  overall  2-D  target-detection  algorithm  based  on  the  coefficient 
change  functional  is  illustrated  in  Fig.  11.  Hidden  within  this  block 

diagram  are  a  number  of  important  decisions.  First,  m  represents  the  local 

s 

mean  of  s(n,m)  under  the  sliding  window  w(n,m).  As  in  the  1-D  case  (see 
Sec.  2.1.5),  in  theory  this  mean  level  should  be  removed,  although  in 
practice  it  was  found  to  make  little  difference  in  the  coefficient  change 
functional.  (However,  the  coefficient  estivates  themselves  may  differ.) 

Thus,  m#  was  not  subtracted  from  the  data,  although  it  was  found  useful  at 
times  as  an  additional  parameter  within  the  detection  functional. 

Another  decision  (hidden  within  the  "direct  estimation"  box)  involves 
the  choice  of  the  shape  and  size  of  the  prediction  mask.  Recall  that  our 
derivations  in  the  previous  sections  were  restricted  to  a  first-quadrant 
quarter-plane  mask.  Since  an  image  has  no  apparent  directionality,  a 
quarter-plane  mask  may  bias  the  coefficient  estimates.  Furthermore,  from  the 
work  of  Ekstrom  and  Woods,7  we  know  chat  a  quarter-plane  mask  (even  infinite 
in  extent)  is  not  sufficient  to  match  sn  arbitrary  power  spectrum,  whereas  a 
nonsymoatric  half-plane  mask  is  sufficient. 
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As  illustrated  in  Fig.  11,  before  applying  a  threshold,  we  smoothed  the 
detection  functional  with  an  FIR  filter.  This  operation  helped  reduce  the 
number  of  false  alarms.  The  threshold  choice  was  based  on  heuristics, 
although  we  are  currently  considering  an  adaptive  threshold  (e.g.,  CFAR  used 
in  radar).  Finally,  to  create  the  target  map,  a  3  x  3  median  filter  was 
applied  to  further  reduce  the  false-alarm  rate. 

The  following  examples  illustrate  these  considerations  and  also 
demonstrate  the  success  of  detecting  targets  of  small  extent. 

2.2.5  Examples 

In  the  following  set  of  examples,  we  have  applied  the  2-D  nonrecursive 
least-squares  algorithm  to  the  RADC  image  data  which  are  composed  of  forests 
and  fields  illustrated  in  Fig.  12.  Anomalous  regions  in  the  photograph  were 
created  by  taking  small  areas  (2  x  2,  3  x  3,  and  A  x  A  pixels  in  extent)  of 
forested  regions  and  transplanting  them  to  the  field  regions  as  depicted  in 
Fig.  13.  The  coefficient-change  technique  can  find  these  anomalous  areas  in 
spite  of  the  use  of  the  very  low-order  model  with  P  ■  Q  ■  1  (i.e,,  three 
independent  coefficients). 

Figure  1A  illustrates  the  coefficient  change  function  Eq.  (45)  (before 
smoothing)  based  on  a  5  x  5  sliding  window4  and  a  2  x  2  first-quadrant 
prediction  mask.  Note  that  C[n,m]  is  largest  at  the  three  targets  and  at  the 
field-field  and  field-tree  transitions.  Figure  15  depicts  c(n,ra]  based  on  a 
2x2  second-quadrant  prediction  mask.  The  average  of  Figs.  1A  and  15  is 
shown  in  Fig.  16  which  appears  to  be  an  improvement  over  either  of  these  two 
functionals,  i.e.,  the  target  and  transition  regions  are  more  clearly 
accentuated.  A  cross  section  of  the  original  image  (through  the 
2x2  target)  and  the  corresponding  cross  section  of  the  averaged  C(n,m)  are 
illustrated  in  Figs.  17(a)  and  (b),  respectively. 

We  might  consider  the  average  in  Fig.  16  to  be  a  rough  approximation  to 
a  functional  based  on  a  nonaymmetric  half-plane  (NSHP)  prediction  mask.7  A 
comparison  of  a  first-quadrant,  aecond-quadrant,  and  a  12-coefficient  NSHP 
prediction  mask  is  illustrated  in  Fig.  18.  C(n,m]  associated  with  this  NSHP 

*A11  examples  in  this  section  are  baaed  on  a  5  x  5  sliding  window. 


Fig.  14,  Detection  functional  baaed  on  f irat-quadraot  meek. 
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Fig.  16.  Average  of  Figs.  14  and  15 


Fig.  17(b).  Cross  section  of  detection  functional  in  Fig.  16. 
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FIRST  QUADRANT 


SECOND  QUADRANT 


NONSYMAftETRiC  HALF-PLANE 

Fig.  18.  Comparison  of  different  prediction  masks. 

mask  performed  about  as  well  as  the  average  depicted  in  Fig.  16  in  detecting 
targets,  although  it  seems  to  exhibit  a  poorer  performance  in  detecting  the 
field-field  transition. 

Figures  19(a)  and  (b)  depict  the  smoothed  C[n,m]  of  Fig.  16  after 
application  of  a  3  x  3  FIR  smoothing  filter.  The  result  of  thresholding  and 
median  filtering  to  create  the  target  map  is  ahovn  in  Fig.  20  superimposed 

on  the  original  image.  All  three  targets  and  the  field-tree  tranaition  are 
clearly  detected. 

These  examples  represent  a  subset  of  a  large  number  of  experiments 
performed  in  2-D  based  on  the  coefficient-change  function.  For  example,  the 
local  mean  level  was  not  used  in  the  experimonts  described  above,  but  was 
found  in  other  examples  tii  enhance  C{n,o]  which  was  reflected  in  reducing  the 
false-alaro  rate  of  the  target  map.  Other  experiments  involved  smoothing  the 
parameters  a(n^,mjJ  with  a  Laplacian  before  computing  Cln^.m^J  (represents  s 
differencing  operation  in  two  directions),  applying  an  exponential  weight  to 
the  windowed  data  before  estimation  (in  order  to  emphasize  certain  regions), 


Fig.  19(a).  Smoothed  detection  functional  for  RADC  image 


Fig.  19(b).  Contour  map  of  detection  functional  in  Fig.  19(a). 
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Fig.  20.  Target  rap  corresponding  to  detection  functional 
and  image  in  Figs.  19  and  13. 
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increasing  and  decreasing  the  window  size,  basing  a  detection  functional  on  a 
geometric  mean  squared  change  (rather  than  an  arithmetic  mean),  incorporating 
a  variance  measurement  into  the  detection  functional,  etc.  One  problem  in 
all  these  studies  involved  the  detection  of  anomalous  areas  in  the  tree  data 
of  Fig.  12.  Preliminary  results  indicate,  however,  that  a  2-D  extension  of 
the  use  of  the  prediction  error  may  have  some  advantages  in  such  detection 
problems  over  the  coefficient-change  function. 
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